Libraries

# SPDS
library(tidyverse)
## Warning: replacing previous import 'vctrs::data_frame' by 'tibble::data_frame'
## when loading 'dplyr'
## ── Attaching packages ───────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ──────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(units)
## udunits system database from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/units/share/udunits
# Data
library(USAboundaries)
library(rnaturalearth)

# Visualization
library(gghighlight)
library(ggrepel)
library(knitr)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(readxl)
library(leaflet)
library(raster) # Raster Data handling
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(tidyverse) # Data Manipulation
library(getlandsat) # keyless Landsat data (2013-2017)
library(sf) # Vector data processing
library(mapview) # Rapid Interactive visualization
## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Question 1: AOI Identification

Palo = read_csv('../data/uscities.csv') %>%
  filter(city == "Palo") %>%
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
  st_transform(5070) %>%
  st_buffer(5000) %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_as_sf()
## Parsed with column specification:
## cols(
##   city = col_character(),
##   city_ascii = col_character(),
##   state_id = col_character(),
##   state_name = col_character(),
##   county_fips = col_double(),
##   county_name = col_character(),
##   county_fips_all = col_character(),
##   county_name_all = col_character(),
##   lat = col_double(),
##   lng = col_double(),
##   population = col_double(),
##   density = col_double(),
##   source = col_character(),
##   military = col_logical(),
##   incorporated = col_logical(),
##   timezone = col_character(),
##   ranking = col_double(),
##   zips = col_character(),
##   id = col_double()
## )

Question 2

Step 2: Read in Metadata

meta = read_csv("../data/palo-flood-scene.csv")
## Parsed with column specification:
## cols(
##   entityId = col_character(),
##   acquisitionDate = col_datetime(format = ""),
##   cloudCover = col_double(),
##   processingLevel = col_character(),
##   path = col_double(),
##   row = col_double(),
##   min_lat = col_double(),
##   min_lon = col_double(),
##   max_lat = col_double(),
##   max_lon = col_double(),
##   download_url = col_character()
## )
files = lsat_scene_files(meta$download_url) %>%
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file)) %>%
  arrange(file) %>%
  pull(file)

Step 3:

st = sapply(files, lsat_image)
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
s = stack(st) %>%
  setNames(c(paste0("band", 1:6)))

The dimensions of the stacked image are 7811 rows by 7681 columns. The Coordinate Reference System used is UTM Zone 15 and the cell resolution is 30 x 30.

Step 4

crop = Palo %>%
  st_transform(crs(s))

The dimensions of the stacked image are 340 rows by 346 columns. The Coordinate Reference System used is UTM Zone 15 and the cell resolution is 30 x 30.

Question 3 – RGB Images

r = crop(s, crop)
par(mfrow = c(2,2))
plotRGB(r, r = 4, g = 3, b = 2, stretch = "lin")
#natural
plotRGB(r, r = 5, g = 4, b = 3, stretch = "lin")
#Traditional Color Infrared
plotRGB(r, r = 5, g = 6, b = 4, stretch = "lin")
#False Color used for Water focus
plotRGB(r, r = 5, g = 7, b = 1, stretch = "lin")
## Warning in .local(x, ...): layer was changed to 6

#uses NIR, SWIR2, Coastal Aerosols used to visualize vegetation in water

Part 2

Applying a color stretch provides greater contrast between the coloration of cells within the raster, and can make it easier to observe desired spatial trends.

Question 4

palette = colorRampPalette(c("blue", "white", "red"))
ndvi = (r$band5 - r$band4) / (r$band5 + r$band4)
plot(ndvi, col = palette(256))

ndwi = (r$band3 - r$band5) / (r$band3 + r$band5)
plot(ndwi, col = palette(256))

mndwi = (r$band3 - r$band6) / (r$band3 + r$band6)
plot(mndwi, col = palette(256))

wri = (r$band3 + r$band4) / (r$band5 + r$band6)
plot(wri, col = palette(256))

swi = (1 / (sqrt(r$band2 - r$band6)))
## Warning in sqrt(getValues(x)): NaNs produced
plot(swi, col = palette(256))

index_stack = stack(ndvi, ndwi, mndwi, wri, swi) %>%
  setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(index_stack, col = palette(256))

Each of the plots help visualize slightly different components of the same area of interest. The NDVI plot shows the vegetated land in red, and non-vegetated land in blue. The NDWI plot shows is a similar plot, but instead shows areas which are water in red and dry-land areas in blue. The MNDWI plot enhances the visualization of the NDWI plot, showing further distinction between water (red) and land (blue). The WRI plot goes even further, normalizing the colors of water and land regions. The final plot completely removes land (shows in white) and leaves us to visualize only the water areas in blue.

Step 2 Thresholding

thresholding1 = function(x){ifelse(x <= 0, 1, NA)}
flood1 = calc(ndvi, thresholding1)
thresholding2 = function(x){ifelse(x >= 0, 1, NA)}
flood2 = calc(ndwi, thresholding2)
thresholding3 = function(x){ifelse(x >= 0, 1, NA)}
flood3 = calc(mndwi, thresholding3)
thresholding4 = function(x){ifelse(x >= 1, 1, NA)}
flood4 = calc(wri, thresholding4)
thresholding5 = function(x){ifelse(x <= 5, 1, NA)}
flood5 = calc(swi, thresholding5)
flood_stack = stack(flood1, flood2, flood3, flood4, flood5) %>%
  setNames(c("NDVI Flood", "NDWI Flood", "MNDWI Flood", "WRI Flood", "SWI Flood"))
plot(flood_stack, col = "blue")

Question 5

set.seed(09102020)
getValues(r) %>%
  dim() %>%
  na.omit()
## [1] 117640      6

The number of extracted values matches the number of cells in the raster, or the dimensions multiplied together. This means that the values were extraceted cell by cell, and that each cell was assigned a single value.

data_r = getValues(r) %>%
  na.omit()
kmean = kmeans(data_r, 12)
fviz_cluster(kmean, geom="point", data = data_r)

thresholding = function(x){ifelse(x <= 0, 1, 0)}
findflood = calc(ndvi, thresholding)
kmeansraster = findflood 
values(kmeansraster) = kmean$cluster
plot(kmeansraster, col = viridis::viridis(12))

com_table = table(values(findflood),values(kmeansraster))
kmeansraster[kmeansraster != which.max(com_table[2,])] = 0
kmeansraster[kmeansraster != 0] = 1
plot(kmeansraster)

thresholdings2 = function(x){ifelse(x >= 0, 1, 0)}
findflood2 = calc(ndwi, thresholdings2)
thresholdings3 = function(x){ifelse(x >= 0, 1, 0)}
findflood3 = calc(mndwi, thresholdings3)
thresholdings4 = function(x){ifelse(x >= 1, 1, 0)}
findflood4 = calc(wri, thresholdings4)
thresholdings5 = function(x){ifelse(x <= 5, 1, 0)}
findflood5 = calc(swi, thresholdings5)
stackfindflood = stack(findflood, findflood2, findflood3, findflood4, findflood5, kmeansraster)
plot(stackfindflood)

sumsff = sum(stackfindflood)
plot(sumsff, col = blues9)

(cellStats(stackfindflood, sum) * res(sumsff)^2) / 1e6
## layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 
##  5.9994  6.4908 10.7451  7.6221 13.6809  7.5780

Extra Credit

AOI = st_point(c(-91.78967, 42.06290)) %>%
  st_sfc(crs = 4326) %>%
  st_as_sf() %>%
  st_transform(st_crs(stackfindflood)) 
raster::extract(sumsff, AOI)
## [1] 2

Two of the rasters identified the point as a flooding area.